Unfallanalyse im Kanton Zürich
Joel Wenger und Christian Wassmer
2022-04-07
Einleitung
Dieser Transferbericht untersucht die Einflussfaktoren von tödlichen Unfällen. Er wurde im Rahmen der Vorlesungen “Lineare Regression” und “Zeitreihenanalyse” der BFH im Herbstsemester 2021/22 erfasst, welche von Dr. Raul Gimeno gehalten wurden.
Datenquellen
Unfalldaten
- https://opendata.swiss/de/dataset/polizeilich-registrierte-verkehrsunfalle-im-kanton-zurich
- https://opendata.swiss/de/dataset/polizeilich-registrierte-verkehrsunfalle-im-kanton-zurich/resource/1b2a47b7-ea16-4264-9bb6-d743cebdbd76
Die Verkehrsunfallstatistik des Kantons Zürich (VUSTA) enthält die Strassenverkehrsunfälle mit Personen- und Sachschäden, die durch die Kantonspolizei Zürich, die Dienstabteilung Verkehr der Stadt Zürich sowie die Stadtpolizei Winterthur registriert wurden. Sie wird einmal jährlich aktualisiert, jeweils gegen Ende des ersten Quartals des Folgejahres.
Zu jedem Strassenverkehrsunfall sind der Unfallort (geokodiert), das Jahr, der Monat, der Wochentag, die Unfallstunde, die Strassenart, der Unfalltyp (ab 1. Juli 2015 inklusive Bagatellunfälle zB. Parkschäden), die Unfallbeteiligung (‘Fussgänger’, ‘Velo (ohne E-Bikes)’, ‘Motorrad’) und die Unfallschwerekategorie verfügbar. Detaillierte Definitionen der Variabeln sind in der Ressource “Minimales Geodatenmodell Strassenverkehrsunfallorte (ASTRA)” dokumentiert.
Die Datei “KTZH_00000718_00001783.csv” wurde am 11.03.2022 heruntergeladen. Am 16.03.2022 wurde eine neue Version veröffentlicht, welche mit dem Script nicht kompatibel ist.
Meteodaten
Der Datensatz umfasst Stundenwerte ab 1992 bis zur letzten aktuellen Stunde aufgeteilt in Jahresdateien. Darin enthalten sind die Stationen Stampfenbachstrasse, Schimmelstrasse und Rosengartenstrasse. Gemessen wird der Luftdruck (p), die Niederschlagsdauer (RainDur), die Globalstrahlung (StrGlo), die Temperatur (T), die relative Luftfeuchtigkeit (Hr), die Windrichtung, die Vektor und Skalar Windgeschwindigkeit. Vor 2018 sind die Skalar Windgeschwindigkeiten aus den 30 Minuten Vektor Daten gerechnet worden. s Die Stundenwerte des laufenden Jahres werden jeweils 30 Minuten nach der vollen Stunde aktualisiert.
Lineare Regression “Unfalldaten”
Daten einlesen
normalizeText <- function(x, length = 9) {
res <- str_trim(x, side = 'both')
res <- str_replace_all(res, ' ', '_')
substr(res, 0, length)
}
# 'basedir' ist eine lokale Variable und verweist auf ein Verzeichnis in meiner Dropbox
accidentsFile <- paste0(basedir, 'KTZH_00000718_00001783.csv')
rawAccident <- read.csv(accidentsFile, header = TRUE, sep = ';', stringsAsFactors = TRUE)
accidentData <- rawAccident %>%
select(1, AccidentType_de, AccidentSeverityCategory_de, RoadType_de, MunicipalityCode, AccidentYear, AccidentMonth, AccidentWeekDay, AccidentHour) %>%
mutate(AccidentType_de = normalizeText(AccidentType_de)) %>%
mutate(AccidentSeverityCategory_de = normalizeText(AccidentSeverityCategory_de, 20)) %>%
mutate(RoadType_de = normalizeText(RoadType_de)) %>%
mutate(YearMonth = paste(AccidentYear, sprintf("%02.f", AccidentMonth), sep = '-'))
# Aus Gründen der Übersichtlichkeit wird die UID des Unfalls ausgeblendet
head(accidentData %>% select(-1))## AccidentType_de AccidentSeverityCategory_de RoadType_de MunicipalityCode
## 1 Schleuder Unfall_mit_Sachschad Hauptstra 247
## 2 Schleuder Unfall_mit_Sachschad Nebenstra 261
## 3 Schleuder Unfall_mit_Leichtver Nebenstra 261
## 4 Schleuder Unfall_mit_Sachschad Autobahn 251
## 5 Schleuder Unfall_mit_Sachschad andere 261
## 6 Überquere Unfall_mit_Leichtver Nebenstra 261
## AccidentYear AccidentMonth AccidentWeekDay AccidentHour YearMonth
## 1 2011 1 aw406 0 2011-01
## 2 2011 1 aw406 0 2011-01
## 3 2011 1 aw406 1 2011-01
## 4 2011 1 aw406 1 2011-01
## 5 2011 1 aw406 2 2011-01
## 6 2011 1 aw406 2 2011-01## Daten einlesenUnfall nach Typ
accidentByType <- table(accidentData$YearMonth, accidentData$AccidentType) %>%
as.data.frame() %>%
pivot_wider(names_from = Var2, values_from = Freq, values_fill = 0) %>%
rename(YearMonth = Var1)
tail(accidentByType, 12)## # A tibble: 12 × 12
## YearMonth Abbiegeun Andere Auffahrun Einbiegeu Frontalko Fussgänge Parkierun
## <fct> <int> <int> <int> <int> <int> <int> <int>
## 1 2020-01 43 18 219 73 13 38 203
## 2 2020-02 64 32 243 79 7 50 205
## 3 2020-03 38 10 147 67 6 23 177
## 4 2020-04 49 12 123 63 7 22 150
## 5 2020-05 70 14 220 86 13 29 189
## 6 2020-06 62 15 249 101 18 48 255
## 7 2020-07 56 22 270 90 23 35 247
## 8 2020-08 56 24 245 73 14 43 257
## 9 2020-09 61 30 261 91 17 43 227
## 10 2020-10 68 23 222 100 22 37 256
## 11 2020-11 69 18 201 70 13 35 204
## 12 2020-12 28 27 195 79 13 38 235
## # … with 4 more variables: Schleuder <int>, Tierunfal <int>, Überholun <int>,
## # Überquere <int>Unfall nach Schwerekategorie
accidentBySeverity <- table(accidentData$YearMonth, accidentData$AccidentSeverityCategory) %>%
as.data.frame() %>%
pivot_wider(names_from = Var2, values_from = Freq, values_fill = 0) %>%
rename(YearMonth = Var1)
tail(accidentBySeverity, 12)## # A tibble: 12 × 5
## YearMonth Unfall_mit_Getö… Unfall_mit_Leic… Unfall_mit_Sach… Unfall_mit_Schw…
## <fct> <int> <int> <int> <int>
## 1 2020-01 0 174 946 26
## 2 2020-02 4 189 977 31
## 3 2020-03 5 138 740 31
## 4 2020-04 4 226 629 39
## 5 2020-05 0 296 910 76
## 6 2020-06 5 304 1033 66
## 7 2020-07 2 345 1056 68
## 8 2020-08 0 299 1048 66
## 9 2020-09 3 323 1044 71
## 10 2020-10 2 223 1109 48
## 11 2020-11 3 198 923 29
## 12 2020-12 2 166 1009 25Unfall nach Strassentyp
accidentByRoadType <- table(accidentData$YearMonth, accidentData$RoadType) %>%
as.data.frame() %>%
pivot_wider(names_from = Var2, values_from = Freq, values_fill = 0) %>%
rename(YearMonth = Var1)
tail(accidentByRoadType, 12)## # A tibble: 12 × 7
## YearMonth Autobahn Autostras Hauptstra Nebenanla Nebenstra andere
## <fct> <int> <int> <int> <int> <int> <int>
## 1 2020-01 173 8 276 2 585 102
## 2 2020-02 182 10 294 1 600 114
## 3 2020-03 121 2 219 0 479 93
## 4 2020-04 72 5 217 3 516 85
## 5 2020-05 132 6 343 2 690 109
## 6 2020-06 150 9 365 2 743 139
## 7 2020-07 185 12 364 4 782 124
## 8 2020-08 191 14 318 1 735 154
## 9 2020-09 197 14 341 4 776 109
## 10 2020-10 162 5 336 5 734 140
## 11 2020-11 148 14 271 2 618 100
## 12 2020-12 142 7 281 0 661 111Unfall nach Stunde
Wird in der Regressionsanalyse nicht weiter verwendet.
accidentByHour <- table(accidentData$YearMonth, accidentData$AccidentHour) %>%
as.data.frame() %>%
pivot_wider(names_from = Var2, values_from = Freq, values_fill = 0) %>%
rename(YearMonth = Var1)
tail(accidentByHour, 12)## # A tibble: 12 × 25
## YearMonth `0` `1` `2` `3` `4` `5` `6` `7` `8` `9` `10`
## <fct> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 2020-01 24 8 8 6 12 20 33 90 63 54 55
## 2 2020-02 20 15 13 10 3 19 54 77 69 50 80
## 3 2020-03 19 10 5 6 4 5 32 45 64 52 61
## 4 2020-04 15 8 4 4 10 5 13 34 32 38 50
## 5 2020-05 18 11 6 7 9 9 34 48 55 57 73
## 6 2020-06 24 17 4 9 8 12 29 65 70 70 94
## 7 2020-07 20 9 7 14 14 20 45 53 57 78 78
## 8 2020-08 21 17 12 10 12 15 37 60 70 61 96
## 9 2020-09 25 11 9 9 11 20 35 73 76 66 81
## 10 2020-10 16 16 12 4 5 13 50 68 72 65 71
## 11 2020-11 20 8 0 7 5 3 55 51 67 58 64
## 12 2020-12 27 20 11 16 8 22 35 50 63 47 74
## # … with 13 more variables: `11` <int>, `12` <int>, `13` <int>, `14` <int>,
## # `15` <int>, `16` <int>, `17` <int>, `18` <int>, `19` <int>, `20` <int>,
## # `21` <int>, `22` <int>, `23` <int>Unfall nach Wochentag
accidentByWeekDay <- table(accidentData$YearMonth, accidentData$AccidentWeekDay) %>%
as.data.frame() %>%
pivot_wider(names_from = Var2, values_from = Freq, values_fill = 0) %>%
rename(YearMonth = Var1) %>%
rename(Mo = aw401, Di = aw402, Mi = aw403, Do = aw404, Fr = aw405, Sa = aw406, So = aw407)
tail(accidentByWeekDay, 12)## # A tibble: 12 × 8
## YearMonth Mo Di Mi Do Fr Sa So
## <fct> <int> <int> <int> <int> <int> <int> <int>
## 1 2020-01 170 158 202 202 180 131 103
## 2 2020-02 190 149 177 189 203 205 88
## 3 2020-03 177 141 137 118 141 104 96
## 4 2020-04 117 116 179 158 128 110 90
## 5 2020-05 172 172 195 193 225 189 136
## 6 2020-06 210 280 247 201 192 183 95
## 7 2020-07 193 213 249 256 244 194 122
## 8 2020-08 228 197 198 185 242 196 167
## 9 2020-09 177 252 266 213 229 193 111
## 10 2020-10 166 181 180 236 270 212 137
## 11 2020-11 192 160 188 145 206 152 110
## 12 2020-12 157 235 223 175 190 137 85Zusammenfassung
Das Datenset erstreckt sich über 120 Monate (Januar 2011 bis Dezember 2022) und enthält rund 150k Unfälle (pro Monat rund 1200).
accidentByMonth <- accidentData %>%
group_by(YearMonth) %>%
summarise(TotalAccidents = n())
tail(accidentByMonth, 12)## # A tibble: 12 × 2
## YearMonth TotalAccidents
## <chr> <int>
## 1 2020-01 1146
## 2 2020-02 1201
## 3 2020-03 914
## 4 2020-04 898
## 5 2020-05 1282
## 6 2020-06 1408
## 7 2020-07 1471
## 8 2020-08 1413
## 9 2020-09 1441
## 10 2020-10 1382
## 11 2020-11 1153
## 12 2020-12 1202NROW(accidentByMonth)## [1] 120mean(accidentByMonth$TotalAccidents)## [1] 1233.092sum(accidentByMonth$TotalAccidents)## [1] 147971min(accidentByMonth$YearMonth)## [1] "2011-01"max(accidentByMonth$YearMonth)## [1] "2020-12"Alle Daten zusammenführen
accidentFreq <- accidentByMonth %>%
# full_join(accidentByHour, by = 'YearMonth') %>%
full_join(accidentByWeekDay, by = 'YearMonth') %>%
full_join(accidentByType, by = 'YearMonth') %>%
full_join(accidentBySeverity, by = 'YearMonth') %>%
full_join(accidentByRoadType, by = 'YearMonth')
# Sanity checks
head(accidentFreq)## # A tibble: 6 × 30
## YearMonth TotalAccidents Mo Di Mi Do Fr Sa So Abbiegeun
## <chr> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 2011-01 1018 141 112 125 244 160 129 107 41
## 2 2011-02 851 127 110 122 119 145 125 103 34
## 3 2011-03 1071 143 192 161 182 177 127 89 46
## 4 2011-04 1159 149 142 184 195 198 186 105 57
## 5 2011-05 1284 217 203 194 165 190 179 136 70
## 6 2011-06 1127 134 174 213 211 149 152 94 52
## # … with 20 more variables: Andere <int>, Auffahrun <int>, Einbiegeu <int>,
## # Frontalko <int>, Fussgänge <int>, Parkierun <int>, Schleuder <int>,
## # Tierunfal <int>, Überholun <int>, Überquere <int>,
## # Unfall_mit_Getöteten <int>, Unfall_mit_Leichtver <int>,
## # Unfall_mit_Sachschad <int>, Unfall_mit_Schwerver <int>, Autobahn <int>,
## # Autostras <int>, Hauptstra <int>, Nebenanla <int>, Nebenstra <int>,
## # andere <int># Still 120 months?
NROW(accidentFreq) == 120## [1] TRUERegressionsanalyse
Hypothese / Untersuchungsgegenstand: Welches sind die Haupteinfluss-Faktoren für “Unfall_mit_Getöteten”?
Mittel: Step-wise Vorwärts- und Rückwärts-Selektion sowie Suche in “beiden Richtungen” mit dem Akaike‘s AIC-Kriterium und dem Schwarz Informationskriterium (BIC).
Vorbereitungen
# Entfernen nicht mehr benötigte Daten
accidents <- accidentFreq %>%
select(-YearMonth, -TotalAccidents, -Unfall_mit_Sachschad, -Unfall_mit_Leichtver, -Unfall_mit_Schwerver)
# Interzeptmodell als Ausgangmodell
nullModel <- lm(Unfall_mit_Getöteten ~ 1, data = accidents)
fullModel <- lm(Unfall_mit_Getöteten ~ ., data = accidents)
N <- NROW(accidents)AIC-Kriterium
forwardModelAIC <- step(nullModel,
direction = "forward",
scope = formula(fullModel), k = 2, trace = 0)
backwardModelAIC <- step(fullModel,
direction = "backward",
scope = formula(fullModel), k = 2, trace = 0)
bothModelAIC <- step(fullModel,
direction = "both",
data = accidents, k = 2, trace = 0)Interpretation
Mit allen drei Verfahren kommen die gleichen Regressoren heraus mit gleichem AIC:
paste(AIC(forwardModelAIC), AIC(backwardModelAIC), AIC(bothModelAIC), sep = ", ")## [1] "448.006256423627, 448.006256423627, 448.006256423627"paste(attr(forwardModelAIC$terms, "term"), sep = ", ")## [1] "Abbiegeun" "Fr" "Fussgänge" "Tierunfal" "Sa" "Autobahn"
## [7] "Frontalko" "So"paste(attr(backwardModelAIC$terms, "term"), sep = ", ")## [1] "Fr" "Sa" "So" "Abbiegeun" "Frontalko" "Fussgänge"
## [7] "Tierunfal" "Autobahn"paste(attr(bothModelAIC$terms, "term"), sep = ", ")## [1] "Fr" "Sa" "So" "Abbiegeun" "Frontalko" "Fussgänge"
## [7] "Tierunfal" "Autobahn"finalModelAIC <- lm(Unfall_mit_Getöteten ~
Fr +
Sa +
So +
Abbiegeun +
Frontalko +
Fussgänge +
Tierunfal +
Autobahn, data = accidents)
summary(finalModelAIC)##
## Call:
## lm(formula = Unfall_mit_Getöteten ~ Fr + Sa + So + Abbiegeun +
## Frontalko + Fussgänge + Tierunfal + Autobahn, data = accidents)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2681 -0.8868 -0.3055 0.9596 3.5694
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.515703 1.054844 2.385 0.018779 *
## Fr -0.015005 0.004392 -3.416 0.000888 ***
## Sa 0.011098 0.005643 1.967 0.051696 .
## So 0.009357 0.006847 1.367 0.174503
## Abbiegeun 0.035209 0.013508 2.607 0.010402 *
## Frontalko 0.057109 0.033439 1.708 0.090458 .
## Fussgänge 0.035559 0.015551 2.287 0.024112 *
## Tierunfal -0.040967 0.021387 -1.916 0.057993 .
## Autobahn -0.019077 0.007067 -2.700 0.008028 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.497 on 111 degrees of freedom
## Multiple R-squared: 0.2517, Adjusted R-squared: 0.1977
## F-statistic: 4.666 on 8 and 111 DF, p-value: 6.103e-05layout(matrix(1:4, 2, 2))
plot(finalModelAIC)
jarque.test(finalModelAIC$residuals)##
## Jarque-Bera Normality Test
##
## data: finalModelAIC$residuals
## JB = 2.5572, p-value = 0.2784
## alternative hypothesis: greaterDie Residuen sind normalverteilt. Das Modell kann verwendet werden, siehe Analyse / Fazit.
BIC-Kriterium
forwardModelBIC <- step(nullModel,
direction = "forward",
scope = formula(fullModel), k = log(N), trace = 0)
backwardModelBIC <- step(fullModel,
direction = "backward",
scope = formula(fullModel), k = log(N), trace = 0)
bothModelBIC <- step(fullModel,
direction = "both",
data = accidents, k = log(N), trace = 0)Interpretation
Mit allen drei Verfahren kommen die gleichen Regressoren heraus mit gleichem BIC:
paste(BIC(forwardModelAIC), BIC(backwardModelAIC), BIC(bothModelAIC), sep = ", ")## [1] "475.881173851447, 475.881173851447, 475.881173851447"paste(attr(forwardModelAIC$terms, "term"), sep = ", ")## [1] "Abbiegeun" "Fr" "Fussgänge" "Tierunfal" "Sa" "Autobahn"
## [7] "Frontalko" "So"paste(attr(backwardModelAIC$terms, "term"), sep = ", ")## [1] "Fr" "Sa" "So" "Abbiegeun" "Frontalko" "Fussgänge"
## [7] "Tierunfal" "Autobahn"paste(attr(bothModelAIC$terms, "term"), sep = ", ")## [1] "Fr" "Sa" "So" "Abbiegeun" "Frontalko" "Fussgänge"
## [7] "Tierunfal" "Autobahn"finalModelBIC <- lm(Unfall_mit_Getöteten ~ Fr + Abbiegeun, data = accidents)
summary(finalModelBIC)##
## Call:
## lm(formula = Unfall_mit_Getöteten ~ Fr + Abbiegeun, data = accidents)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3389 -0.9213 -0.1741 0.7858 4.7462
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.002103 0.773064 2.590 0.010822 *
## Fr -0.010850 0.003845 -2.822 0.005611 **
## Abbiegeun 0.049589 0.013203 3.756 0.000271 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.584 on 117 degrees of freedom
## Multiple R-squared: 0.1167, Adjusted R-squared: 0.1016
## F-statistic: 7.732 on 2 and 117 DF, p-value: 0.0007015layout(matrix(1:4, 2, 2))
plot(finalModelBIC)
jarque.test(finalModelBIC$residuals)##
## Jarque-Bera Normality Test
##
## data: finalModelBIC$residuals
## JB = 6.0343, p-value = 0.04894
## alternative hypothesis: greaterDie Residuen sind nicht normalverteilt. Das Model wird nicht weiter verfolgt.
Analyse und Fazit
Uns ist bewusst, dass mit der vorliegenden Methode keine belastbaren, geschweige denn kausalen Zusammenhänge hergeleitet werden können. Zumal die Kennzahl R2 mit rund 0.25 darauf hinweist, dass das gewählte Regressionsmodell den Regressanden “Unfall_mit_Getöteten” nur ungenügend schwach erklären kann. Dennoch scheint sich zu bewahrheiten, dass die Wochentage Fr (***), Sa (.) und So (nicht signifikant), sowie Abbiegeunfälle (*), Frontalkollisionen (.), Unfälle mit Fussgängern (*) und auf Autobahnen (**) den stärksten Einfluss auf Unfälle mit Todesopfern haben. Dies entspricht unserer Ansicht nach durchaus den Erwartungen. Einzig der Einflussfaktor “Unfälle mit Tieren” (.) kommt etwas unerwartet.
(In Klammern Signifikanzcodes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1)
Das Modell, welches mittels BIC gefunden wurde, wurde verworfen: Die Residuen sind nicht normalverteilt.
Alle Aussagen basieren auf einem Signifikanzniveau von 95%.
Interessant wäre vermutlich die Einflussfaktoren in den Gruppen einzeln zu analysieren, also Wochentage, Unfallart, Strassentyp und ggf. auch Stunde; als auch mit anderen Schwerekategorien zu untersuchen.
Zeitreihenanalyse “Unfall- und Meteodaten”
Daten einlesen
Unfalldaten
Es werden dieselben Daten wie bereits für das Regressionsmodell verwendet.
# function to change column names (remove "Accident", whitespaces and language designator)
process_headers <- function(Headers) {
NewHeaders <- NULL
for (i in 1:length(Headers)) {
tmp <- str_replace(Headers[i], '_', '')
tmp <- str_replace(tmp, 'Accident', '')
NewHeaders[i] <- str_replace(tmp, 'en', '')
}
return(NewHeaders)
}
# function to extract calendar day from available data
process_day <- function(Month, Weekday) {
Day <- 1
for (i in 2:length(Weekday)) {
# increment day counter if weekday variable changes
if (Weekday[i] == Weekday[i - 1]) {
Day[i] <- Day[i - 1]
} else {
Day[i] = Day[i - 1] + 1
}
# reset if month variable changes
if (Month[i] != Month[i - 1]) {
Day[i] = 1
}
}
return(Day)
}
# load raw data into workspace
AccidentFile <- paste0(basedir, 'KTZH_00000718_00001783.csv')
RawAccident <- read.csv(AccidentFile, header = TRUE, sep = ';', stringsAsFactors = TRUE, colClasses = c("AccidentUID" = "character"))
# process raw data to extract information needed for this project
Accident.df <- RawAccident %>%
select(6, 11, 12, 13, 14, 19, 20, 21, 23, 24, 25, 30, 35) %>%
drop_na %>%
rename_with(process_headers) %>%
mutate(WeekDay = as.numeric(substr(WeekDay, 5, 5))) %>%
mutate(Day = process_day(Month, WeekDay)) %>%
mutate(InvolvingPedestrian = ifelse(InvolvingPedestrian == 'true', 1, 0)) %>%
mutate(InvolvingBicycle = ifelse(InvolvingBicycle == 'true', 1, 0)) %>%
mutate(InvolvingMotorcycle = ifelse(InvolvingMotorcycle == 'true', 1, 0)) %>%
mutate(IsSevere = ifelse(grepl('severe', SeverityCategory) | grepl('fatalities', SeverityCategory), 1, 0)) %>%
mutate(IsDeadly = ifelse(grepl('fatalities', SeverityCategory), 1, 0))Meteodaten
Die Meteodaten werden direkt von data.stadt-zuerich.ch heruntergeladen und in einen dataframe gespeichert. Wurde die Datei bereits einmal ausgeführt, so ist MeteoDataZurich.Rda bereits verfügbar und wird lokal geladen.
# check if data exists in directory
if (!file.exists(paste0(basedir, 'MeteoDataZurich.Rda'))) {
# create empty dataframe
MeteoData.df <- data.frame()
# get values for 2011 - 2021
for (i in 2011:2021) {
url <- paste("https://data.stadt-zuerich.ch/dataset/ugz_meteodaten_stundenmittelwerte/download/ugz_ogd_meteo_h1_", as.character(i), ".csv", sep = '')
data <- read.csv(url, encoding = 'UTF-8')
df <- data.frame(data) %>%
mutate(Year = as.numeric(substr(X.U.FEFF.Datum, 1, 4))) %>%
mutate(Month = as.numeric(substr(X.U.FEFF.Datum, 6, 7))) %>%
mutate(Day = as.numeric(substr(X.U.FEFF.Datum, 9, 10))) %>%
mutate(Hour = as.numeric(substr(X.U.FEFF.Datum, 12, 13))) %>%
pivot_wider(names_from = Parameter, values_from = Wert) %>%
select(Year, Month, Day, Hour, T, RainDur) %>%
group_by(Year, Month, Day, Hour) %>%
summarise(Temperature = mean(T, na.rm = TRUE), RainDuration = mean(RainDur, na.rm = TRUE)) %>%
ungroup %>%
select(Year, Month, Day, Temperature, RainDuration) %>%
group_by(Year, Month, Day) %>%
summarise(Temperature = mean(Temperature, na.rm = TRUE), RainDuration = sum(RainDuration, na.rm = TRUE))
MeteoData.df <- rbind(MeteoData.df, df)
}
# save data to Rda file
save(MeteoData.df, file = paste0(basedir, 'MeteoDataZurich.Rda'))
} else {
load(paste0(basedir, 'MeteoDataZurich.Rda'))
}Daten zusammenführen
# add meteo information to dataframe
Accident.df <- merge(Accident.df, MeteoData.df, by = c('Year', 'Month', 'Day'))Daten bearbeiten und gruppieren
Tagesdaten
# group data into daily sets
AccidentDaily.df <- Accident.df %>%
select(Year, Month, Day, InvolvingPedestrian, InvolvingBicycle, InvolvingMotorcycle, IsSevere, IsDeadly, Temperature, RainDuration) %>%
group_by(Year, Month, Day) %>%
mutate(Total = n()) %>%
summarise(Total = mean(Total),
InvolvingPedestrian = sum(InvolvingPedestrian),
InvolvingBicycle = sum(InvolvingBicycle),
InvolvingMotorcycle = sum(InvolvingMotorcycle),
IsSevere = sum(IsSevere),
IsDeadly = sum(IsDeadly),
Temperature = mean(Temperature),
RainDuration = mean(RainDuration)) %>%
mutate(Date = as.Date(paste(as.character(Year), '-', as.character(Month), '-', as.character(Day), sep = ''), format = '%Y-%m-%d')) %>%
select(Date, everything()) %>%
ungroupMonatsdaten
# further group data into monthly sets
AccidentMonthly.df <- AccidentDaily.df %>%
group_by(Year, Month) %>%
summarise(Total = sum(Total),
InvolvingPedestrian = sum(InvolvingPedestrian),
InvolvingBicycle = sum(InvolvingBicycle),
InvolvingMotorcycle = sum(InvolvingMotorcycle),
IsSevere = sum(IsSevere),
IsDeadly = sum(IsDeadly),
Temperature = mean(Temperature),
RainDuration = sum(RainDuration)) %>%
mutate(Date = as.Date(paste(as.character(Year), '-', as.character(Month), '-01', sep = ''), format = '%Y-%m-%d')) %>%
select(Date, everything()) %>%
ungroupGrafische Datenexploration
Histogramme
# accidents vs. hour
hist(Accident.df$Hour)
hist(Accident.df$Hour[Accident.df$InvolvingPedestrian == TRUE])
hist(Accident.df$Hour[Accident.df$InvolvingBicycle == TRUE])
hist(Accident.df$Hour[Accident.df$InvolvingMotorcycle == TRUE])
# accidents vs. weekday
hist(Accident.df$WeekDay)
hist(Accident.df$WeekDay[Accident.df$InvolvingPedestrian == TRUE])
hist(Accident.df$WeekDay[Accident.df$InvolvingBicycle == TRUE])
hist(Accident.df$WeekDay[Accident.df$InvolvingMotorcycle == TRUE])
# accidents vs. month
hist(Accident.df$Month)
hist(Accident.df$Month[Accident.df$InvolvingPedestrian == TRUE])
hist(Accident.df$Month[Accident.df$InvolvingBicycle == TRUE])
hist(Accident.df$Month[Accident.df$InvolvingMotorcycle == TRUE])
Korrelationen
# daily data
CorrDataDaily.df <- AccidentDaily.df %>%
select(Total, InvolvingPedestrian, InvolvingBicycle, InvolvingMotorcycle, IsSevere, IsDeadly, Temperature)
CorrMatrixDaily <- cor(CorrDataDaily.df)
ggcorrplot(CorrMatrixDaily, type = "lower", lab = TRUE)
#ggpairs(CorrDataDaily.df) # attention -> takes a very long time to visualize (large dataset)
# monthly data
CorrDataMonthly.df <- AccidentMonthly.df %>%
select(Total, InvolvingPedestrian, InvolvingBicycle, InvolvingMotorcycle, IsSevere, IsDeadly, Temperature, RainDuration)
CorrMatrixMonthly <- cor(CorrDataMonthly.df)
ggcorrplot(CorrMatrixMonthly, type = "lower", lab = TRUE)
ggpairs(CorrDataMonthly.df)
Zeitreihen
Zeitreihen plotten
# select some data to make ready to plot with ggplot
AccidentMonthly.plot.df <- AccidentMonthly.df %>%
select(Date, InvolvingPedestrian, InvolvingBicycle, InvolvingMotorcycle) %>%
pivot_longer(!Date, names_to = 'Involved', values_to = 'Values')
# plot data using ggplot
ggplot(AccidentMonthly.df, aes(x = Date, y = Total)) +
geom_line() +
ylab('Number of Accidents') +
ggtitle('Total Number of Monthly Accidents') +
theme_light()
ggplot(AccidentMonthly.plot.df, aes(x = Date, y = Values, color = Involved)) +
geom_line() +
ylab('Number of Accidents') +
ggtitle('Monthly Number of Accidents involving pedestrians, bicycles and motorcycles') +
theme_light()
Zeitreihenanalyse
# create timeseries from monthly data
InvolvingBicycleMonthly.ts <- ts(AccidentMonthly.df$InvolvingBicycle, start = c(AccidentMonthly.df$Year[1], AccidentMonthly.df$Month[1]), frequency = 12)
InvolvingMotorcycleMonthly.ts <- ts(AccidentMonthly.df$InvolvingMotorcycle, start = c(AccidentMonthly.df$Year[1], AccidentMonthly.df$Month[1]), frequency = 12)
# decompose time series objects
InvolvingBicycleMonthly.ts.m <- decompose(InvolvingBicycleMonthly.ts, type = 'multiplicative')
plot(InvolvingBicycleMonthly.ts.m)
InvolvingMotorcycleMonthly.ts.a <- decompose(InvolvingMotorcycleMonthly.ts, type = 'additive')
plot(InvolvingMotorcycleMonthly.ts.a)
Analyse und Fazit
Grundsätzlich entspricht die Zeitreihenanalyse den Erwartungen, welche man aufgrund der Saisonalität von Unfällen mit Fahrrädern und Motorrädern haben könnte: In den wärmeren Monaten, insbesonderen an warmen Tagen (starke Korrelation zwischen Temperatur und Anzahl Velo- und Motorradunfällen), sind die Unfallzahlen mit Velos und Motorrädern stark erhöht gegenüber den Wintermonaten. Die Zahl der Velounfälle scheint auch mit einem gewissen Trend zu steigen.
Zudem scheint eine negative Korrelation mit Regentagen und Unfällen mit Velos und Motorrädern zu bestehen, was ebenfalls Sinn machen würde.
Im Allgemeinen sieht man nichts Aussergewöhnliches, als das was man aufgrund der Jahreszeiten sowieso vermuten würde.